### code largely based on SI material for Borcard D., Gillet F. & Legendre P. 2018.
## Numerical Ecology with R, 2nd edition. Springer International Publishing AG.  
## and course material generously provided by A. Buttler (EPFL Numerical Ecology Class)

# 
# comment -----------------------------------------------------------------
# blabla
# this is a comment
# ctrl + shift + R 

# data exploration --------------------------------------------------------




# sub section -------------------------------------------------------------



# help --------------------------------------------------------------------
help(vegan)
?vegan
??vegan


# working directory -------------------------------------------------------
getwd()
setwd(dir = "C:/Users/hpeter.INTRANET/Desktop/multivariate stats R/2024/scripts/")
setwd(choose.dir())


# Install, load packages and data -----------------------------------------
#install.packages("vegan")
library(vegan)
library(ade4)

# load data ---------------------------------------------------------------
load("Doubs.RData")
data (doubs)

# read data from file -----------------------------------------------------
ex1<-read.csv("example.csv", head=TRUE, sep=",", row.names = 1)
ex2<-read.csv(file.choose())


# define data
a=12
b=34
c=a+b


dat1<-c(2,3,5)
dat1*2



# Exploring a dataset ------------------------------------------------
spe                       # Display the whole data frame in the console 
X<-spe[c(2,5,7), 8:12]            # Display only 5 lines and 10 columns
head(spe)                 # Display only the first 6 lines
tail(spe)                 # Display only the last 6 rows
nrow(spe)                 # Number of rows (sites)
ncol(spe)                 # Number of columns (species)
dim(spe)                  # Dimensions of the data frame (rows, columns)
B=colnames(spe)             # Column labels (descriptors = species)
rownames(spe)             # Row labels (objects = sites)
summary(spe)              # Descriptive statistics for columns
str(spe)                  # Structure of the dataset


# subsetting a dataframe --------------------------------------------------
ablette<-fishtraits[fishtraits$EnglishName=="Bleak",]
body_length2<-fishtraits$BodyLength
body_length=fishtraits[,6]
body_length_ablette=fishtraits[fishtraits$EnglishName=="Bleak",6]
body_length_ablette


# descriptive statistics using iris dataset --------------------------------------------------
data(iris)
head(iris)
summary(iris[iris$Species=="versicolor",])
mean(iris$Sepal.Length)
sd(iris$Sepal.Length)
median(iris$Petal.Width)
var(iris$Sepal.Length[iris$Species=="virginica"])

attach(iris)
mean(Sepal.Length)
hist(iris$Sepal.Length)
detach(iris)
detach()

#install.packages("pastecs")
library(pastecs)
?stat.desc
stat.desc(iris)
stat.desc(iris)[11,2]

#install.packages("psych")
library(psych)

?describeBy

attach(iris)
describeBy(iris, group=Species)
describeBy(iris, group=Species)[3]
detach(iris)



# distribution of abundance (doubs) ---------------------------------------

# Minimum and maximum of abundance values in the whole data set
range(spe)

# Minimum and maximum value for each species
apply(spe, 2, range)  #apply a function (range) to a dataset (spe) over columns (MARGIN=2) or rows (MARGIN=1)
?apply
apply(spe,1,range)



# Count the cases for each abundance class
ab <- table(unlist(spe))
ab

# Barplot of the distribution, all species confounded
barplot(ab, 
        las = 1,
        xlab = "Abundance class",
        ylab = "Frequency",
        col = gray(5 : 0 / 5),
        horiz=F,
)

# Number of absences
sum(spe == 0)

# Proportion of zeros in the community data set
sum(spe == 0) / (nrow(spe) * ncol(spe))

# matrix visualization
library(gplots)
heatmap.2(as.matrix(spe), dendrogram = "none", trace="none")



# Map of the locations of the sites ===============================

# Create an empty frame (proportional axes 1:1, with titles)
# Geographic coordinates x and y from the spa data frame
plot(spa, 
     asp = 1, 
     type = "n", 
     main = "Site Locations", 
     xlab = "x coordinate (km)", 
     ylab = "y coordinate (km)"
)
# Add a blue line connecting the sites along the Doubs River
lines(spa, col = "blue", lwd=3)
# Add the site labels
text(spa, row.names(spa), cex = 1, col = "red")
# Add text blocks 
text(68, 20, "Upstream", cex = 1.2, col = "red")
text(15, 35, "Downstream", cex = 1.2, col = "red")




# Maps of some fish species =======================================

# Divide the plot window into 4 frames, 2 per row
par(mfrow = c(2,2))
# Plot four species
plot(spa, 
     asp = 1, 
     cex.axis = 0.8, 
     col = "brown", 
     cex = spe$Satr, 
     main = "Brown trout", 
     xlab = "x coordinate (km)", 
     ylab = "y coordinate (km)"
)
lines(spa, col = "light blue", lwd=3)

plot(spa, 
     asp = 1, 
     cex.axis = 0.8, 
     col = "brown", 
     cex = spe$Thth, 
     main = "Grayling", 
     xlab = "x coordinate (km)", 
     ylab = "y coordinate (km)"
)
lines(spa, col = "light blue", lwd=3)

plot(spa, 
     asp = 1, 
     cex.axis = 0.8, 
     col = "brown", 
     cex = spe$Baba, 
     main = "Barbel", 
     xlab = "x coordinate (km)", 
     ylab = "y coordinate (km)"
)
lines(spa, col = "light blue", lwd=3)

plot(spa, 
     asp = 1, 
     cex.axis = 0.8, 
     col = "brown", 
     cex = spe$Abbr, 
     main = "Common bream", 
     xlab = "x coordinate (km)", 
     ylab = "y coordinate (km)"
)
lines(spa, col = "light blue", lwd=3)


# Compare species: number of occurrences ==========================

# Compute the number of sites where each species is present
# To sum by columns, the second argument of apply(), MARGIN, 
# is set to 2
spe.pres <- apply(spe > 0, 2, sum)
spe.pres
# Sort the results in increasing order
sort(spe.pres)

# Compute percentage frequencies
spe.relf <- 100 * spe.pres/nrow(spe)

# Round the sorted output to 1 digit
round(sort(spe.relf), 1)

# Plot histograms
# Divide the window horizontally
par(mfrow = c(1,2)) #  ?par
hist(spe.pres, 
     main = "Species Occurrences", 
     right = FALSE, 
     las = 1, 
     xlab = "Number of occurrences", 
     ylab = "Number of species", 
     breaks = seq(0, 30, by = 5),
     col = "bisque"
)
hist(spe.relf, 
     main = "Species Relative Frequencies", 
     right = FALSE, 
     las = 1,
     xlab = "Frequency of occurrences (%)", 
     ylab = "Number of species",
     breaks = seq(0, 100, by = 10),
     col = "bisque"
)




# Compare sites: species richness =================================

# Compute the number of species at each site
# To sum by rows, the second argument of apply(), MARGIN, is set 
# to 1
sit.pres <- apply(spe > 0, 1, sum)
# Sort the results in increasing order
sort(sit.pres)

par(mfrow = c(1, 2))
# Plot species richness vs. position of the sites along the river
plot(sit.pres,type = "s",
     las = 1, 
     col = "blue",
     main = "Species Richness vs. \n Upstream-Downstream Gradient",
     xlab = "Site numbers", 
     ylab = "Species richness"
)
text(sit.pres, row.names(spe), cex = .8, col = "red")

# Use geographic coordinates to plot a bubble map
plot(spa, 
     asp = 1, 
     main = "Map of Species Richness", 
     pch = 21, 
     col = "white", 
     bg = "brown", 
     cex = 5 * sit.pres / max(sit.pres), 
     xlab = "x coordinate (km)", 
     ylab = "y coordinate (km)"
)
lines(spa, col = "light blue", lwd=3)


vignette(package="vegan")
vignette("diversity-vegan")

# HOMEWORK: calculate and display different alpha-diversity indices using the vegan R package!


# Environmental Data ======================================================

par(mfrow = c(2, 2))
plot(env$dfs, env$ele,
     xlab = "Distance from the source (km)", 
     ylab = "Elevation (m)", pch=16,
     col = "red", main = "Elevation"
)
plot(env$dfs, env$dis, 
     type = "l", 
     xlab = "Distance from the source (km)", 
     ylab = "Discharge (m3/s)", 
     col = "blue", 
     main = "Discharge",
     lwd=2
)
plot(env$dfs, env$nit, 
     type = "o", 
     xlab = "Distance from the source (km)", 
     ylab = "Nitrate (mg/L)", 
     col = "brown", 
     main = "Nitrate"
)


plot(env$dfs, env$oxy, 
     type = "b", 
     xlab = "Distance from the source (km)", 
     ylab = "Oxygen (mg/L)", 
     col = "green3", 
     main = "Oxygen", 
    cex=env$nit
)


#relationship between Nitrate and Oxygen
par(mfrow=c(1,1))
plot(env$oxy~env$nit,  pch=16, xlab="Nitrate (mg/L)", ylab="Oxygen (mg/L)")
mod<-lm(env$oxy~env$nit)
summary(mod)
abline(mod, lwd=2, lty=3)



# radar plots
library(fmsb)
summary(env)

max_min <- data.frame(
  dfs = c(450, 0), ele = c(1000, 150), slo = c(50, 0),
  dis = c(70, 0), pH = c(9, 6), har = c(120, 20))
rownames(max_min) <- c("Max", "Min")
env2 <- rbind(max_min, env[,1:6])
env2

sample4<-env2[c(1,2,6),]

radarchart(sample4)  #radar chart single site


par(mar = c(1, 2, 2, 2)) #radar chart multiple sites
radarchart(env2[1:6,])
legend(x="bottomright", legend = rownames(env2[3:6,]), horiz = TRUE,
       bty = "n", pch = 20, col=1:4, title="sites")




# Scatter plots for all pairs of environmental variables ==========
pairs(env)
pairs(env[1:5])

?pairs



panel.hist <- function(x, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(usr[1:2], 0, 1.5) )
  h <- hist(x, plot = FALSE)
  breaks <- h$breaks; nB <- length(breaks)
  y <- h$counts; y <- y/max(y)
  rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}

panel.cor <- function(x, y){
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- round(cor(x, y), digits=2)
  txt <- paste0("R = ", r)
  cex.cor <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = cex.cor * r)
}


pairs(env[1:5], 
      panel = panel.smooth, 
      diag.panel = panel.hist,
      upper.panel=panel.cor,
      cex.labels=3,
      main = "Bivariate Plots with Histograms and Smooth Curves"
)

panel.cor2 <- function(x, y){
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- round(cor(x, y), digits=2)
  txt <- paste0("R = ", r)
  text(0.5, 0.5, txt, cex = 2)
}


pairs(env[1:5], 
      panel = panel.smooth, 
      diag.panel = panel.hist,
      upper.panel=panel.cor2,
      cex.labels=3,
      main = "Bivariate Plots with Histograms and Smooth Curves"
)